#include<iostream>
using namespace std;
#include<cmath>
#include<fstream>
#include"library.h"

double E;

void derivative(int n, double x, double *y, double *dydx){
	dydx[0] = y[1];
	dydx[1] = - E* y[0];
}

double psib(double x){

	E=x;

	int N=10000,n=2;
	double a=0, ya=0, deriv_a=1, b=1;
	int i,j;
	double h = (b-a)/N;
	double *y = new double [n];
	double *ynext = new double [n];
	y[0]=ya; y[1] = deriv_a;
	for(i=0;i<N;i++){
		rk4(a+i*h, h, n, y, ynext, derivative);
		for(j=0;j<n;j++) y[j] = ynext[j];
	}
	return y[0];
}

int main(){

	double eig = bisection(psib, 50, 100, 1E-9);
	cout<<"E/pi/pi="<<eig/M_PI/M_PI<<endl;

	return 0;
}
